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NATHBMATICAL  PROBUBMS  IN  CREATINO  A  DIGITAL  GRID  ANALXZER 


N.A.  Aleksldze 

A  basic  obstacle  to  achieving  great  accuracy  In  solving  boundary 
problems  for  equations  In  mathematical  physics  by  the  grid  method  Is 
the  ccanparatlvely  small  memory  of  present  electronic  di^i^ltal  compu¬ 
ters.  As  a  result.  It  Is  necessary  to  divide  the  grid  Into  parts  dur¬ 
ing  solution  of  these  problems,  and  to- return  to  the  external  storage 
upon  each  Iteration,  or  to  use  Schwartz's  alternating  method.  As  was 

shown  In  note  [2],  the  latter  method  Is  preferable  when  the  grid  re^* 

* 

glon  Is  divided  Into  two  parts.  If  the  number  of  subsets  of  mesh 
points  to  be  added  Is  Increased,  the  rate  of  convergence  of  Schwartz's 
Iteration  process  decreases  rapidly,  and  this  method  too  becomes 
rather  Inefficient.  We  should  note  that  when  the  region  Is  divided  ! 
Into  parts,  programming  difficulties  appear  as  well.  In  view  of  what  | 
we  have  said.  It  Is  a  very  attractive  Idea  to  construct  a  specialized  I 
digital  computer  whose  working  storage  uses  either  a  magnetic  drum  ' 
[14]  or  closed  magnetic  tape  [16],  which  could  be  used  for  problems 
Involving  Laplace,  Poisson,  wave,  or  diffusion  eq^atlons.  The  mathe- 

I 

matloal  prerequisite  for  this  Idea  is  slaqpllolty  of  carrying  through  I 
the  iteration  process,  both  with  the  very  simple  and  with  the  nine- 
point  grid  approximation  of  the  Laplaolan,  which  enters  Into  all  the 
equations  mentioned  above. 

The  aim  of  ths/prement  article  is  to  investigate  the  deslrablULty 
from  tte  vieiqpolat  of  maehine  time,  of  using  various  iterat'lon  pro- 
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cesses  in  specialized  digital  computers.  As  we  shall  show  below,  the 
parallel  processing  proposed  earlier  for  all  columns  (rows),  while 
substantially  ^nlarglng  the  apparatus /yields  almost  no  advantage 
with  the  most  esqpedient  iteration  process,  the  superrelaxation  proc-  * 


In  order  to  investigate  the  esqpedlency  of  using  the  various  iter¬ 
ation  processes,  it  is  necessary  to  know  the  as3rnq;>tptlo  relationships 
between  the  local  and  Integral  errors  for  the  different  iteration 
processes.  These  relationships  are  derived  In  §1. 

In  §2,  we  examine  the  various  iteration  processes  from  the  view¬ 
point  of  convenience  of  the  calculations  set  up  and  machine  tine. 

In  S3i  we  discuss  one  grid  approximation  for  the  Laplacian  opera¬ 
tor,  and  one  effective  iteration  process  for  solution  of  the  first 
boundary  problem  is  proposed. 

In  $4,  we  derive  relationships  between  the  congmter  digit  capa9|^ 

Ity  and  the  size  of  the  working  storage  that  are  optimum  from  the 

viewpoint  of  a  difference  solution  to  the  boundary  problems. 

§1.  Asymptotic  Relationships  Between  Local  and  Integral  Errors  in 
Solving  Boundary  Problems 


1.  We  know  that  the  difference  equations  obtained  by  substituting 
difference  relationships  for  elliptical  differential  operators  sui^e,  as 
a  rule,  solved  by  an  Iteration  method.  Only  residues  can  be  used  to 
find  the  end  of  the  Iteration  process  and  the  Integp^al  error  in  the 
approximate  solution  to  the  difference  equations.  Nonetheless,  the 
literature  contains  no  discussions  of  probldns  of  the  asyoqptotio  rela¬ 
tionship  between  the  residue  and  the  Integral  error  in  various  itera¬ 
tion  processes.  Oershgorin  [8]  derived  an  accurate  estimate  of  the  . 
error  in  tezms  of  the  residue  for  the  solution  of  difference  equatims 
by  the  majorant  sMthod.  Per  the  slapleet  grid  approxlBatio^  of  tiM  La- 


plaolan  differential  operator,  the  Gershgorin  estimate  takes  the  fol¬ 


lowing  form 


While  for  the  nine -point  approximation 

I  I  - 


(1.1) 


(1.2) 


Where  le  the  maxlTOiTn  residue,  1  Is  the  radius  of  a  circle  Inter¬ 
secting  the  given  region  as  a  whole,  h  Is  the  mesh  Increment.  Zt  Is 
clear  from  (1.1),  for  example,  that  where  l/h  >  60,  the  residue  should 
be  three  orders  of  magnitude  smaller  than  the  anticipated  accuracy 

—if 

e  -  10  ,  and  this  cannot  be  ensured  if  the  machine  contains  fewer 
than  (k  +  3)  decimal  digits,  to  say  nothing  of  the  fact  that  this 
leads  to  protraction  of  the  Iteration  process.  In  this  paragraph,  we. 
shall  derive  asynqptotlc  estimates  for  the  various  Iteration  processM, 
allowing  fairly  large  residues  In  order  to  x)rovlde  the  required  accu^ 


racy. 


2.  We  write  the  finite -difference  equivalent  of  the  Dlrlchlet 


problem  for  the  Laplace  equation  In  matrix  form 


(1.3) 


where  t  ^  0  only  at  the  boundary  mesh  points.  For  an  error  of  e  ■  u  — 
-  u,  the  equation  Ac  ■  A(u  —  u)  ■  or  e  ■  A“^o^  j  will  hold, 
where  o.  .  Is  the  residue  at  the  (l,j)  mesh  point  and  u  Is  an  approx- 
Imate  solution  to  Bq.  (1.3).  We  write  A  In  the  following  form  [12] i 


■ 


whore  D  ■  B  -i-  rA,  B  Is  the  unit  ,  matrix,  r  Is  a  relaxation  fadtor  so 
chosen  that  the  spectral  nonuD  of  operator  B  .«•  rP  i^U  be  a  mlnfma. 


{ 

» 
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Since  r  may  be  so  chosen  that  all  characteristic  values  of  operator  D 
will  have  modu^  less  than  unlty^  when  n  ea^  0«  and  for  e«  we 


obtain 


1*1 


(1.4) 


We  let  ph  and  qh  be  the  height  and  width  of  the  smallest  rectangle 
completely  enclosing  the  given  region.  For  the  slnqplest  approximation 
of  the  two-dimensional  Laplaclem  operator^  r,  and  the  corresponding 
estimate  e  will  take  on  the  loLlowlng  values: 

For  Richardson's  Iteratloa  process,  r—  —  ,  Omi—  -- (y*  +^"*) 

4  4 


(1.5) 


For  Llebman's  Iteration  process,  Dm  “(P“*+  f"*)J  and 


(1.6)  k 


JVTWT- 


(1.7) 


For  a  second-order  Richardson  process  [13]«  £  must  be  taken  eus 
the  sum  of  the  optimum  values  of  the  two  parameters  that  determine 
the  rate  of  convergence  for  this  process.  Then  r  <  1,  2hn — 


For  a  superrelaxatlon  Iteration  process  (Llebman's  method  extrap* 
olated  [13])  Dmi—  YT  ;  and 


|.|T  . 


(1.8) 


Estimates  (1.3) and  (1.6)  permit  larger  residue  values  for  ensiu*- 
Ing  the  required  accuracy  than  does  the  Oershgorln  estimate  (1.1)  for 
any  region,  while  the  linear  dimensions.  In  \anit8  of  h,  appear  In  ee- 
tlmates  (1.7)  and  (1.8)  In  first  degree,  idiloh  substantlAlly  reduoes 
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the  number  of  iterations. 

3*  The  ohlu'aoteristlc  values  of  the  nine -point  approximation  M  t<^ 
the  Laplaolan  operator  that  have  the  largest  and  smallest  moduli  are,i 
for  the  Richardson  iteration  process,  equal  to,  respectively  [51 «  I 

I  A/.— 

while  we  obtain 

I  ^27  a .  t 


for  the  optimum  value  of  r. 

For  this  value  of  r,  the  spectral  norm  of  the  operator  | |E  +  rM| 
equals  »  <uid  for  the  estimate  e  we  obtain 


(1.9) 


A  Richardson  iteration  process  of  second  order  is  carried  out 
for  the  operator  M  in  accordance  with  the  following  formula 

If  r^  and  rg  are  calculated  from  the  equations 

i—ft  iM-'.+a/'r^-c^ 

then  we  obtain  »  which  leads-  to  the  following  esti¬ 

mate  for  e 

4.  Estimates  (1.5) -(1.10),  Just  as  in  the  case  of  the  Oershgorin 
estimates  (1.1)  and  (1.2),  cannot  be  used  if  drift  operatmm  az>e  used 
for  the  boundary  values  in  the  region.  Let  \as  examine  the  following 
difference  operators* 


*  0  ■> 


1 


where 


C*  I  *•+»•  >+tJ  2  (l -;-&(+ /4.j) 


I  ^  I  «  \'  **-f* 


*+^+r/+f 

•f<4»»/»^]  ■>■  . ^ - fTTT^ -  (*!+»./« 

^  +•»!+»»  i+«)j — «fci» 


M“0.  ±‘»  ^<“0.  ±»5  r*+i*>o; 


6^^^  j  Is  the  distance  from  mesh  point  (l,J)  to  the  boundary  r  In  the 
positive  direction  with  respect  to  the  Ox  axis  (to  point  A),  divided 
by  h  (Fig.  l);  j  the  value  of  function  u  at  point  A  on  bound- 

ary  T;  Is  the  distance  from  mesh  point  (1,J)  to  the  bounda^ 

r  In  the  direction  of  the  bisector  of  the  angle  formed  by  lines  x 
-  Xj^  €md  y  ■  y^  (to  point  B),  divided  by  is  the  value 

of  function  u  at  point  B  on  the  boundary  r,  etc. 
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Th6  following  two  lemmas  were  pi-oved  in  [!i]: 

Lemma  1. »  l: 

4 

ferenoe  equation 


Lemma  1,  *  if  la  the  maxlmvun  modulus  of  the  solution  to  a  dlf- 


DKus^muj,  .  fkjmak, 
Ftoy-o,  a<ri^iA. 

i.  J'-®' 


tlon 

M  ®Wi 

ir|r-iO, 

ruiUih, 

• 

and  satisfies  the  Inequality 

. ..  . . . V. .  _ -  _ 

(1.11) 

Lemma  2.  If  M..  Is 

the  maximum  modulus  of  the 

solution  to  the 

difference  equation 

BY 

i  F4./—0, 

F,-o, 

r*/SaA, 

.  o<r(i./<aA, 

then  H.  Is  the  maximum  modulus  to 

the  solution  of 

the  difference  equa- 

tlon 

8,.^ 

r^/Kaa, 

o<ffc/<aA, 

l/Ir-O 

and  satisfies  the  Inequality 

Mia 


(1.12) 


(i,+i)(i+i)(i+«  r - „i^ 

where  r^  j  le  the  minimum  dl'etanoe  from  mesh  point  (l«j)  to  the  bound¬ 
ary  of  the  region, ^•“'^f^“*M*»/+J*  ^"^C>jy»*«^^i-iBatlnatea  (1.11) 
and  (1.12)  are  aoourate  inaamuch  as  they  are  also  valid  for.  the  aaya^- 


1  '•nr 

totlo  estimates.  Substituting  estimates  (1. 3)>(1. 10)  Into  (1.11)  auid 
(1.12)  In  plaqe  of  vie  obtain  an  estimate  for  the  error  where  the 
operators  D  and  or  D  and  Mg  are  employed  simultaneous ly^  for  the 
various  Iteration  processes. 

3.  The  boundary  problem  for  the  rectangle  p  «•  30«  q  ■  32  was 
solved  by  the  following  two  methods 

(A) 

where  r  =  0.452  is  the  optimum  value  of  the  releucatlon  factor  for  the 
given  rectangle.  In  both  cases,  the  Iteration  processes  were  contin¬ 
ued  until  the  residues  became  smaller  than  10“^.  The  first  method  re¬ 
quired  986  Iterations,  the  second  method  86.  The  maximum  Integral 

error  In  the  first  case  was  reached  at  the  center  of  the  rectangle,  , 
-4 

and  equaled  2> 10  ,  and  In  the  second  case  It  was  found  at  the  mesh 

point  (1  ■  15«  j  ■  8),  and  equaled  2.2>10'‘^.  The  results  are  In  good 
agreement  with  estimates  (1.6)  and  (1.7)*  If  the  Oershgorln  estimate 
(1)  were  used  to  obtain  an  Integral  error  of  2.2*10'’^,  the  residue 
would  have  to  have  been  reduced  to  less  than  6. 2*  10 which  would 
have  required  Iterations  of  type  (A)  and 

mi6  Iterations  of  type  (B).  Consequently,  In  ccmiparlng 
the  convergences  of  various  Iteration  processes.  In  general  the  resi¬ 
dues  should  not  be  reduced  to  the  same  size  for  all  the  Iteration 
processes,  since  the  same  yields  different  Integral  errors  for  the 
various  Iteration  processes.  We  should  also  note  that  In  Iteration 
process  (B),  the  approximation  already  Is  not  aperiodic;  thus,  for 
example.  If  the  boundary  problem  Is  solved  for  u  s  1,  and  ■  0  Is 
tedcen  as  a  first  anproxinatlon,  while  the  Iteration  process  (B)  is 
progressing,  values ..u  >  1  may  an^ear  at  several  mesh  points,  which 


-8 


must  be  taken  Into  account  in  solving  problems  on  a  fixed-point  c(»qpu- 
ter.  \ 

S2.  Iteration  l^ocesses  Connected  with  a  Difference  Solution  to  the 
Flrsi  Bounoary  Problem 

» 

In  oaloulatlng  naohlne  tlne'^  we  shall  assume  that  the  working 
storage  uses  a  magnetic  drum  with  a  capacity  of  23,000  numbers,  and  a 
speed  of  50  rps. 

The  author  of  [14]  deals  basically  with  the  Richardson  Iteration 
process 

I  - — - 

. _ _  _ 

In  order  to  examine  the  desirability  of  this  process,  we  shall  discuss 
the  two-dimensional  Olrlchlet  problem  for  a  square  with  side  150  h 
(22,300  mesh  points).  We  shall  assume  that  the  Initial  residues  are 

For  the  solution  It  Is  necessary  that  the  Integral  errors  be  less 
than  10 As  was  shown  in  $1,  In  order  for  the  Richardson  Iteration 

process  to  provide  an  Integral  error  of  10*^  In  the  squeo'e  \mder  con-  ' 

I 

slderatlon.  It  Is  necessary  ttot  tlM  looail  errors,  1.  e. ,  the  residues, 
be  brought  to  a  magnitude  of  — jjjs — •o.®*  ‘®  -*  The  spectral  norm  H  of  : 
operator  Ti  equals  K  ■  1  -  it^/2*  130^,  while  the  number  of  Iterations 
required  Is  N» 

vldually,  then  one  Iteration  will  require  I30  cycles,  and  the  solution 

I 

will  take  about  86  hours.  If  all  of  the  columns  have  processed  In  j 

parallel,  about  130  sequential -type  single-digit  accumulators  will  be  | 

, 

required,  and  the  calculation  time  will  be  about  40  minutes.  The 
Richardson  method  Is  clearly  best  suited  to  magnetic  tape,  since  ttie 
writing  of  a  new  value  into  an  old  position  creates  actional  diffl* 


la  [0.8 » 10* 


>io}.5oo.  If  each  column  Is  processed  Indl- 


oultlM.  On  magnetic  tape«  however «  It  le  difficult  to  set  up  such  a 
laxge  number  (130)  of  tracks.  11ius«  we  shall  not  discuss  magnetic 
tapes  In  the  future. 

For  the  Llebman  Iteration  process 

the  spectral  norm  of  operator  A  for  the  problem  under  consideration 

-0 

will  be  Jf— [i— ■•/a*!!©*]*,  the  local  errors  must  be  reduced  to  1.7*10 
and  we  obtaln^^  *  jlii(i^7«io~*3  j  ^  fox*  the  number  of  Iterations.  With  ' 
sequential  processing  of  the  columns «  this  will  take  about  42  hours. 

With  parallel  processing  of  the  columns ,  a  Llebman  Iteration  process 
cannot  be  carried  out«. since  the  (k  +  l)th  approximation  of  the  func¬ 
tions  at  mesh  point  (l^J-l)  is  still  imknown  when  the  (k  +l)th  approx-  , 
Imatlon  Is  being  computed  for  the  function  at  mesh  point  (l«j).  In  ' 
this  oase«  as  a  result,  only  the  following  Iteration  process  Is  possi¬ 
ble 


L  +«&  )• 

4  ■  '  ■  . 


(2.1) 


In  order  to  find  the  spectral  norm  for  operator  A  with  Iteration 
process  (2.1),  we  examine  the  following  Iteration  process 

(2.2) 

which  Is  the  same  as  (2.1)  when  o  -  1/4. 

Let  us  seek  the  characteristic  functions  In  the  form 

(2.3) 


iiW— lia 

15©  *50 


We  write  (2.2)  as  follows: 

Substltutii«  (2*3)  into  (2.4)«  we  obtain 


(2.4)' 


»  M  • 


jnj  .  JtfJ  Tifl  ..  .  nr/  TiSI 

Ale  tia - SID  -^A’sin - sia  — '-4- 

150  *5°  '  ‘S®  *5° 


+.[*»■  d.4:£t:l .1.  ri 

L  150  ISO  150 

+*■•1.  J!±  H,  «'(/-») 

^  ISO  150  •  250  i$o 

150  J  V.  150/  ‘SO 

+ .  I jp«  .i.  i:»±L>  +  .i.  Jii 

_ L _ J  ‘s® 


■ID  J1£L+ 


|-4Jt<aia^  siB 
i$o 


ISO 


We  find  the  arbitrary  constant  R  from  the  condition  A  <■  R  ;  we 
then  obtain 

a  «— I  +4  «— a 


where 


/»— eoa  ; 
_  .15®  _ _*S®_ 


Where  a  ■  l/4«  the  largest  characteristic  value  corresponds  to 
mesh  point  r  >  s  -  !«  and  equals 


Iteration  Process  (2. 1}  will  In  this  case  converge  4/3  times 
faster  than  Richardson's  Iteration  process,  and  2/3  times  more  slowly 
than  Llebmcm's  Iteration  process.  With  parallel  processing  of  all 
columns,  and  with  Iteration  Process  (2. 1)  carried  out,  the  machine 
time  will  be  about  25  minutes. 

The  Richardson  and  Llebman  Iteration  processes  have  an  Inqportant 
advantage  over  other  Iteration  processes  In  that  they  do  not  require 
multiplication  or  division  (division  by  4  Is  accomplished  by  a  shift 
during  write-in). 

The  computer  must  be  capable  of  multiplication,  however.  In  order 
to  process  boundary  mesh  points  for  an  arbitrary  regl<m,  and  In  order 
to  solve  boundary  problems  In  the  case  of  selfad Joint  differential 

-  11  -  ■  ■ 
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I  '*rr 
erators. 

We  shall  now  discuss  the  superrelaxation  Iteration  process 

for  suoh  a  machine.  The  optimum  value  of  the  relaxation  factor  a  is 
found  from  the  equations 


i$o 


-4«+i-o 


and  Is  approximately  equal  to 


X 

300 


As  close  as  this  value  is  to  one-half,  systematic  use  In  (2.5)  of 
the  value  a  -  l/2  leads  to  divergence  of  the  Iteration  process.  The 
spectral  norm  for  the  operator  |  |E  -  ci^p^Al  |  equals 

75 

If  a  superrelaxatlon  Iteration  process  Is  used  to  solve  the 
problem  \uider  consideration,  the  residues  must  be  reduced  to  values  of 
0. 8* 10 and  the  number  of  Iterations  will  be 


la  4 


-350. 


A  solution  will  take  about  I8  minutes  when  the  columns  are  pro¬ 
cessed  sequentially. 

A  superrelaxatlon  Iteration  process  with  parallel  processing  of 
the  columns  Is  Impossible.  Consequently,  we  shall  examine  the  Itera¬ 
tion  process 

with  the  optimum  paraiaeter  a. 

An  inoraaae  in  a  from  1/4  to  1/2  will  decrease  the  oharaoterlstio 
value  A  (r  •  1,  §  m  1)  at  mesh  point  (r  -  1,  a  «  1),  but  thp  modulus 

■  -  12  - 


NiBllSMWBfc 


■  , 

of  the  complex  characteristic  value  at  the  mesh  point  r  »  1,  s  -  149 
will  rapidly  begin  to  increase 

i  (f»i;  V  a*/f4-4« 

where 


t  i 


At  a  ■  l/3i  the  modulus  of  the  characteristic  value  A  (r  ■  1,  s 
■  149)  is  already  the  spectral  norm  of  the  operator 

The  optimum  a  is  calculated  from  the  following  equetion: 

!  el,rfK  * 

Squaring  both  sides  of  this  equation,  we  obtain  a  third -degree 


equation  in  a 


The  discriminant  D  of  this  equation  equals 
■  "  «+>«)*  .  4«+j6)  . _ » 


a4(ad4+4<J)— (4+i6)* 


9(a«4+4i»  ] 


and  can  be  shown  to  be  negative,  which  leeuls  to  inconveniences  in  the 
use  of  Cardan's  formulas.  In  contrast  to  the  superrelaxation  method, 
on  the  other  hand,  an  accurate  optimum  value  for  a  will  in  this  case 
yield  an  insignificant  Increue  in  the  rate  of  convergence.  Since 
where  o  ■  3/10  the  spectral  norm  for  the  operator  | |S  -  0.3A| |  -  7 


will  equal 


and  the  iteration  process  converges  9/7  times  more  rapidly' than  when 
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I  ^ 


a  -  1/4. 

Thus,  whe^  iteration  Process  (2.6)  Is  executed  with  an  optimum 

relaxation  parameter,  the  machine  time  will  be  about  17  minutes. 

< 

A  second-order  Rlchai^dson  Iteration  process  has  a  convergence  of 
about  the  same  order  as  a  superrelaxation  process,  but  Is  Inconvenient 
to  carry  out,  since  It  requires  the  two  preceding  approximations  to 
find  a  new  approximation  and,  consequently,  twice  as  large  a  memory. 

Finally,  we  can  say  that  the  most  expedient  Iteration  process 
for  a  specialized  digital  computer  is  superrelaxation  with  sequential 
processing  of  columns.  If  the  computer  cannot  execute  the  multiplica¬ 
tion  operation,  then  It  Is  best  to  execute  Iteration  Process  (2.3)  us¬ 
ing  a  >  1/4  and  a  >  l/2  alternately. 

S3.  One  Grid  Approximation  for  the  Laplaclatn  Differential  Operator 
Let  us  consider  the  following  finite -difference  equivalent  for 
the  Laplaclan  differential  operator  L 

4  «<*>' 

Let  us  examine  the  various  iteration  processes  using  the  operator  s. 

Let  us  assume  that  the  solution  Is  sought  In  a  rectangle  measur¬ 
ing  ph  and  qh.  We  use  the  following  formula  for  the  Richardson  Itera¬ 
tion  process: 

The  characteristic  functions  of  the  operator  s  will  then  be 

'.if#.*, 

t  1 

V«  now  oaloulate  the  corresponding  characteristic  values  ; 

.is  +ri.  2iiiz±+ 

P  t  P  1- 

+d:^S±:iL  d.  ■■’.'TO, 


t.  ^ 


-4*ln-il..in  JiiL 


M  (  COS  Itl  -cos  —  —  I  )  sin  iCl  sin  fii; 

?  /  P  q  . 


The  spectral  norm  of  this  operator  will  be 

w  909  -y-COB  —  x  J»-8^fa«*(^-»+^-»). 

One  of  the  minimum  characteristic  values  is  obtained  for  the 
point  r  a  s  a  Ij  and  equals 

The  optimum  value  for  a  is  obtained  from  the  condition  1  -f  as^  a 
a  .(i  4-  os*  and  equals  1/4.  The  spectral  norm  for  operator  B  —  (l/4)s 
for  the  Richardson  iteration  process  is 

/-*-4  <>'*+«■*>• 

For  the  Liebman  iteration  process,  the  working  formula  is  written 


we  look  for  the  characteristic  functions  in  the  form 


f  f 


and  we  then  will  have 


*  r 

;+i|Wsis  -fiiiZiL'j. 

.  .  '  f  J  A  I  t 

o 

Selecting  the  coefficient  A  as  follows:  A  a  s,  we  obtain 


- - . __  A  .  f 
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t 


For  the  characteristic  'u_ut  .■  s,  we  have 

W-— 4  «  U-{-4  Ct—l  -  0, 


where 


.  ISf 

faaeOI  -  •CO* 

•  P 


%s 
1  * 


If  we  again  set  a  ■  1/4  (ordinary  Llebman  Iteration  process),  we  ob¬ 
tain 


—'COS*  — 

.  ^  1  _ 

and  the  spectral  norm  s  for  the  operator  E—  (l/4)s  will  eq^al 


If  a  Is  so  chosen  that  A  has  multiple  roots  (Llebman  method  extrapo¬ 
lated),  we  obtain  for  a  the  condition 

4  «*<*««— 4  «+i«o, 


whence 

I  w  - - 

•-  -^~v  r*+?’V » 

»-4  y  r*+r*  • 

The  second -03?der  Richardson  Iteration  process  Is  based  upon  the 
analogy  between  the  solution  to  the  Dlrlchlet  problem  and  the  equation 
d^u/(d  +  2)  +  ^/dt  —  Lu  ■  0,  whose  finite -difference  equivalent  may 
be  written,  for  carrying  out  Integration,  as 


where  r^  and  r2  are  relaxation  factors  whose  optimum  values  are  calcu¬ 
lated  from  the  equations  [13] 

i+r«»*+ri— 


These  values  of  the  relaxation  factors  yield  for  the  spectral  iw>nB  of 
the  oomspoadlnt  operator  a  value 


Similar  problems  for  the  simplest  grid  approximation  to  the  La- 
plaoian  operator  D 

were  discussed  in  [13]. 

C<miparison  of  the  spectral  norms  for  the  operators  shows  that 
the  Richardson  and  Liebman  iteration  processes  converge  twice  as  fast 
for  the  operator  b  as  they  do  for  the  operator  0,  while  the  extrapo¬ 
lated  Liebman  method  and  the  Richardson  iteration  process  of  second 
order  converge  more  rapidly  by  a  factor  of  It  is  not  difficult  to  \ 
show  that  the  operator  s  yields  two  Independent  systems  of  equations 

r-o,>;  i-b,  y; 

••a  iH"> » ^ 

where  k  is  an  Integer.  We  have  used  this  property  of  the  operator  s  to  | 
reduce  the  number  of  arithmetic  operations  in  the  iteration  process. 

In  order  to  check  the  results  obtained,  the  Dirlchlet  problem  was 
solved  for  the  Laplace  equation 

in  a  rectangle  with  sides  30h  and  32h,  using  the  following  methods: 

I.  Liebman  iteration  process  with  operator  D 

I 

II.  Liebman  iteration  process  with  operator  s 

i 

III.  Superrelaxation  iteration  process  with  operator  D  (the  opt!-; 
mum  relaxation  factor  equaled  0.452). 

IV.  Superrelaxation  iteration  process  with  operator  s  (optimum 
relaxation  factor  equaled  0.461). 

V.  Superrelaxation  iteration  pz>ooess  with  operator  s  for  syetem 


(a)  alone,  then  a  single  application  of  operator  D  to  system  (B) 
alone,  and  then  method  III. 

The  value  of  the  function  on  the  boundary  fCs)  varied  within 
wide  limits  (four  problems  were  solved  In  all),  but  the  number  of  It- 
eratlons  for  all  of  the  problems  nonetheless  remained  constant  to 
within  The  calculations  were  halted  when  the  local  errors  became 
less  than  2“^^.  Ve  should  note,  however,  that  as  $1  Implies,  such  a 
local  error  results  In  a  different  Integral  error  for  the  various 
Integration  processes.  Since  the  operator  £  Is  somewhat  less  accurate 
than  the  operator  D,  the  solution  obtained  with  method  XI  was  cor¬ 
rected  by  means  of  method  I,  and  method  IV  by  method  III.  The  corres¬ 
ponding  numbers  of  Iterations  required  for  such  a  correction  are  given 
In  parentheses  In  Table  1. 

TABIE  1 


I 

«  1 

IK 

IV 

<IIICM  KTCptWlA 

WCM  apReMCT* 

IS«7 

»3»(9) 

«3 

71(a) 

7i(>) 

Tm.  onepuiHa 

i6a>io* 

•S'lo* 

i63>io* 

pO'IO* 

1)  Number  of  Iterations;  2)  number  of 
arithmetic  operations. 

The  number  of  Instructions  necessary  to  obtain  a  new  approxima¬ 
tion  for  each  mesh  point  Is  12  with  methods  I  and  II,  and  14  for  me¬ 
thods  III-IV-V.  The  fact  that  the  difference  in  the  nuxnber  of  Instruc¬ 
tions  Is  so  small  Is  erqplalned  by  the  fact  that  the  residue  was  calcu¬ 
lated  In  all  oases  In  order  to  determine  the  point  at  which  the  Itera¬ 
tion  process  was  to  halt. 

GOTparison  of  these  data  with  theoretical  values  Indicated  good 
ngmimentj  alttwiigh  with  methods  UI-IV,  a  sooswhat  smaller  mniber  of 

itermtieni  were  pNMtteted*  The  eeldulaHens  were  eemed  out  <ni  the 
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BBSM  computer. 

54.  Electronic  Dlpcltal  Computer  Parameters  that  nphim™  from  the 
•Tleifpolnt^of  §oundw4roblem  Difference  'Solutions - 

I 

In  order  to  design  a  specialized  digital  computer  to  solve  bound¬ 
ary  problems  [l4j  16]  we  may  derive  the  most  favorable  relationship 
between  the  size  of  the  computer  operating  storage  and  the  number  of 
digits  in  machine  numbers.  In  this  case,  naturally,  the  indicators  of 
the  differential  properties  for  that  class  of  functions  to  which  the 
solution  for  the  problem  under  consideration  belongs  will  enter  into 
these  relationships. 

Let  us  examine  the  general  n-dimenslonal  second -order  differen¬ 
tial  equation  of  elliptical  tyx>e 

.  **  ■  '*•  14.1) 


dll  * 

+*•(*1. (4f|. 

With  the  very  general  boundary  conditions 


V  rf/ 

- _ 


0» 


(4.2) 


and  let  us  write  s<»ne  finite-difference  equivalent  for  the  problem 
(4.  1)-(4.2) 


SW), 


(4.3) 


where  ^  is  the  number  of  adjacent  points,  in  the  sense  of  the  differ¬ 
ence  operator  p^  used;  N  is  the  number  of  mesh  points  in  the  closed 
n-dlmenslonal  region  D  for  idiich  the  problem  (4.1),  (4.2)  is  stated. 
For  ^  0  and  Ag  ^  0,  the  points  at  which  the  lines  generating-  the 
grid  region  intersect  the  boundary  of  the  regi<m  are  also  r^carded  as 
mesh  i>oints.  * 
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I  : 

As  a  rule.  System  (‘t.  j)  lo  .■  .'a.  .01.  method.  It  can 

be  shovm  that,  when  the  conditions  g  <  0,  <  0,  ai.d  not  simultan¬ 

eously  identically  equal  to  0  are  fulfilled ^  the  simple  iteration 
process  for  System  (4.3)  converges.  Let  us  write  the  general  Iteration 
formula  for  (4.3) 

^  ^  Giiy  b2»  -  C„h-cA  u<i>  I4. 

(/-»,»,•••,  AO, 

. _ _ 


Where  r.  and  a  .(s  >1^  2,  ...«  1)  are  the  peuTameters  of  the  Itera- 

J  BpJ  •• 

tlon  formula^  and  are  so  selected  that  the  convergence  will  be  as 

rapid  as  possible;  1  Is  the  number  of  different  approximations  of  the 

« 

function  that  are  used  to  reinforce  convergence.  It  is  assumed  that 


the  (k  +  l)th  approximation  of  the  function  is  written  at  location 
(k-  1),  and  that  all  quantities  entering  into  (4.4),  with  the  excep¬ 
tion  of  u.  ,,  are  computed  once  at  the  beginning  of  the  Iteration 

X  j  J 

process  and  then  remain  stored  In  the  memory.  It  Is  not  difficult  to 
see  that  iteration  Formula  (4.4)  will  then  require  (q  +  e  +  3)  storage 
locatloxis  for  each  mesh  point  (each  number  occupies  one  location).  If 
we  let  be  the  number  of  locations  In  the  operating  storage  of  the 

con^ter,  and  N2  be  the  number  of  Instructions  In  the  Iteration  pro¬ 
gram  [4],  wo  obtain  (Nj^  -  N2)/(q  +1  +  3)  for  the  number  of  mesh 


points.  With  no  loss  of  generality,  we  can  consider  D  to  be  a  unit  n- 

(Nt—Nt 

‘  I 

assume  that  |u-u|  * M  )“*'*  *  the  well- 

known  Banach  theorem  on  convergence  for  the  method  of  successive  approx* 
imatlons,  the  residue  In  Iteration  Prooess  (4.4)  will  reach  a 


value[( 


LfMtjJ 

Lj 

in  order  to  provide  an  aooura^ 
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I  M-rr 


e,  >rtiere  | |p| |  Is  the  spectral  norm  for  the  corresponding  difference 
operator. 

We  must  also  take  Into  acoount  rounding  errors  of  two  types:  a) 
the  error  e'  due  to  rounding  In  the  course  of  Iteration  Process  (4.4); 
b)  the  error  e”  due  to  rounding  In  calculating  the  coefficients  for 
Eq.  (4.3). 


The  mathematical  expectation  for  the  mean-square  error  e'  In  the 
case  of  the  Dlrlchlet  problem  for  the  Laplace  equation  [1] 
after  the  kth  Iteration  due  to  all  of  the  preceding  roundings  will 
first  go  up  as  In  k^  and  will  then  approach  asymptotically 


sin  * '  Fr(^  this  It  follows  that  for  the  residues  to  be 

(JJ 

tJ — 1  j  ,  rounding  must  be  carried  out*,  in 


wins 

reduced . to  a  value 


the  right-hand  digit  of  the  number 


9+H"}  / 


/J 


It  Is  easy  to  verify  that  e"  ■  u  -  fl  satisfy  the  following  dif¬ 
ference  equation  (u  and  Gt  are  solutions  to  (4.3)  for  the  accurate  and 
approximate  values  of  the  coefficients,  respectively) 


where  ((6)  Is  the  rounding  error  In  6.  The  mathematical  expectation 
for  4(6)  equals  zero  but,  nonetheless,  the  probability  estimate  may 
yield  a  value  that  Is  too  low  for  e"  for  the  case  In  which 
g(x^,  ...,  x^)  >  const,  and  f(x2»  ••.«  x^)  const.  The  error  e"  does 
not  depend  upon  the  method  used  to  obtain  a  solution  to  (4.3). 

We  shall  now  show  that  for  the  two-dimensional  Poisson  equation, 
and  for  the  Ll'ebman  and  Richardson  Iteration  processes. 


?H^L, _ 

(the  syittbol  Rmca  indicates  that  I  o<sm  4*^<**)* 

*,  .  I  _ 

-  21  - 


(4.5) 


t 


Let  U8  find  Qreen's  difi’e.en;- 


ij.  following 


(  8  ^ 

manner:  R'  *  ^  for  fixed  Indices  r  and  s  la  a  solution  to  the  follow¬ 


ing  difference  equations: 


*‘'^>1,-0. 


The  solution  of  the  difference  equations 


(4.5') 


can  be  represented  In  terms  of  the  Qreen  differenoe  function  R 


(P*8) 


the  following  manner: 


•5^-EE 

For  In  [5l  we  have  the  estimate 

R‘Cra{x+8)* 


(4.6) 


(4. 7) 


where  a  Is  the  side  of  a  square  completely  enclosed  by  the  given  re¬ 
gion  (the  estimate  given  In  [11]  Is  used  In  deriving  (4. 7))>  Substi¬ 
tuting  (4.7)  Into  (4.6)  and  assuming  that  the  number  of  mesh  points  is 


of  the  order  of  h  ,  we  obtain 


k*lnh 


(4.8) 


where  B  is  some  constant  that  does  not  depend  upon  h.  For  the  Liebnan 

O 

and  Richardson  iteration  processes,  ||p||^l-  oh  and,  therefore. 


**  ^  ^ 

,  »~W  **i»*^‘. 


(4.9) 


where  o^  also  does  not  depend  upon  h.  We  obtain  (4.5)  from  (4.8)  and 
(4.9).  ' 


^  8a  - 


I 


Finally^  we  obtain  for  the  optimum  number  of  digits  in  the  man¬ 
tissa  of  the  number  for  fixed  N 

- - - •, 


(where  2  is  the  radix  of  the  number  system  used  in  the  oonqputer.  We 
assume  that  by  transforming  Eqs.  (4.1) -(4. 2),  the  solution  of  the  dlf. 
ference  equations  will  be  |u|  <  1  (for  example,  for  the  Dlrlchlet 
problem,  we  can  use  the  principle  of  the  maximum  to  do  this).  This 
Imposes  definite  limitations  upon  the  coefficients  of  Eqs.  (4. 1) , 
(4.2). 

Let  us  examine  acme  concrete  cases:  -  Ng  -  10,CXX),  y  ■  2,  a 

superrelaxation  Iteration  process  Is  used,  a^  s  1,  b^^  5  0,  (k  -  1, 

2,  . . . ,  n) , 

yao,  A^SAtSo,\M\at,  o-comt-— (SU  » 
then  we  obtain  frtmi  (4.10) 


mml 

•-a 

•-3 

i 

36 

a4  • 

* 

»3 

1 

3» 

In  practice,  evidently,  we  most  frequently  encounter  problems  In 
which  the  boundary  function  f  ^  H(P,  A,  X),  where  P  +  X  <  2.  For  the 
Dlrlchlet  problem,  with  the  Laplace  equation,  the  best  estimates  are 
contained  In  [6,  7].  In  particular.  It  Is  shown  In  [6]  that  the  error 
is  0(h^^)  for  the  simplest  approximation  to  the  Laplace  operator, 
with  p  -I-  X  <  2.  It  follows  from  [9]  that  asynqptotloally  (h  O) ;  the 
order  of  these  estimates  cannot  be  Improved.  It  has  been  shown  In 
[15]  that  If  a  discontinuity  of  the  boundary  function  occurs  at  a  mMh 


point,  the  error  will  be  .  e  auLoii  processes 

from  the  point  of  view  of  the  number  of  iterations  and  simplicity  of 
the  apparatus  needed  for  the  calculation  scheme  is  the  superrelaxa- 
tlon  iteration  process.  It  follows  from  what  we  have  said  that  when 
|n|  <  1,  for  a  specialized  digital  computer  of  the  types  discussed  in 
[14,  l6],  the  relationship  between  <y^p^  and  N  will  be  (Ng  ■  0,  since 
these  oonputers  use  a  schematic  program) 

L  kIbs  Ina  ~  In  a  J 

In  this  case,  it  is  assumed  that  the  algorithm,  given  in  [3]  is  used 
for  the  Poisson  equation. 

If  no  provision  for  rounding  is  made  in  the  computer,  then  the 
mathematical  expectation  for  the  erroi*  e  'due  to  aii^^ai'diiig  of  the 
right-hand  digits  will  satisfy  Eq.  (4.  5)>  we  derive  the  estimate  e 
in  a  manner  similar  to  that  used  for  e",  we  obtain  for  such  a  compu¬ 
ter 


u 

to/'"'-"*')! 

.  \9H+3). 

1 

_ _ _  _  ...  Iny 

In  conclusion,  we  wish  to  thank  Ye.  A.  Volkov  for  his  helpful  comments 
and  advice. 
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4  Inequality  (1.4)  also  follows  at  once  from  the  familiar 
Banach  convergence  theorem  for  the  method  of  successive 
approximations j  if  A  is  regarded  as  an  operator  that 
transforms  the  space  into  Itself.  In  the  case  of  a 
space  m^«  [tA[|  »  but  within  this  space  the  condition 
||A|J  <  1  is  sufficient  but  not  necessary  for  convergence 
of  the  successive  approximations. 

5  Ml  and  VU  are  the  drift  operators  for  the  function 
boundary'^values  according  to  Mlkeladze  [10],  where  one 
or  two  neighboring  points,  respectively,  are  missing  at 
the  boundary  point.  In  this  notation,  however,  there  is 
a  possibility  that  the  boundary  values  will  drift  not 
only  along  the  lines  forming  the  grid  region,  but  also 
along  lines  parallel  to  the  bisectrices  of  the  coordinate 
angles  (for  lp|  =  |ql  =  1).  They  are  essentially  non- 
symmetrlcal  difference  approximations  of  the  Laplacian 
differential  operator  and  have  a  local  error  of  the  order 
of  l(h3}  under  certain  conditions. 

7  In  Lemmas  1  and  2,  the  angles  at  which  r^  j  <  2h  are  re¬ 

garded  as  boundaries.  Such  "expansion"  of  the  boundary- 
point  zone  is  carried  out  because  the  following  premise 
can  be  demonstrated:  If  the  boundary  of  a  region  is  a 
closed  Lyapunov  curve,  then  there  exists  an  hQ  >  0  such 

that  for  any  h  ^  hQ  at  any  mesh  point  where  r^  j  <  2h  Mj^ 
or  Mg  nay  be  written  in  at  least  one  notation  so  that 

^1+p,  J+q  ^  ”l  *1+P,  J+q  ^^*52  *i+r,j+o 

Mn.  The  proofs  of  Lesmas  1  and  2  do  not  rely  on  this 
assumption  and  if  more  rigid  requirements  are  imposed  on 
the  boundary  of  the  region  and  the  boundary  belt  is 
"narrowed"  to  r^  <  h,  then  6  ^  1  ^  6,  which  improves  the 

evaluations  (l.ll)  and  (1.12). 


* 
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